# 2022.3.24
# edit in 2023.7.26

# Load libraries
library(Seurat)
library(tidyverse)
library(tidyseurat)
library(tidySingleCellExperiment)
library(ggpubr)
library(ggrepel)
library(scMerge)

# load seurat objects
sobj <- read_rds('CRC-I/data/zy2020_tumor10x.rds')

sg_sobj <- read_rds('CRC-I/data/seekgene/crc-starsolo-annotated.rds')

sg_sobj <- sg_sobj |>
  mutate(genotype = case_match(orig.ident, 'crc221125' ~ 'IT', .default = 'TT'))

sg_tam <- sg_sobj |>
  filter(zhang2020_fine == 'TAM-C1QC')

zh_tam <- sobj |>
  filter(str_detect(Sub_Cluster, 'TAM-C1QC'))

gene_merged <- tibble(gene = rownames(sg_tam)) |>
  filter(gene %in% rownames(zh_tam)) |>
  pull(gene)

# visualize CD45 expression as an ref to see if expression scale is the same between 2 data sets.
mrg_tam <- merge(zh_tam, sg_tam)

# keep genes appeared in both datasets
mrg_tam <- subset(mrg_tam, features=gene_merged)

# try scMerge ---------
data("segList", package = "scMerge")

system.time(scm_sce <- mrg_tam %>%
  as.SingleCellExperiment() %>%
  scMerge::scMerge(ctl = segList$human$human_scSEG,
                   assay_name = 'scMerge',
                   hvg_exprs = "logcounts",
                   kmeansK = c(1, 1),
                   batch_name = 'Platform',
                   verbose = TRUE))

scm_sce@assays@data@listData[["scMerge"]][1:5,1:5]

scm_sobj <- scm_sce@assays@data@listData[["scMerge"]] |>
  CreateSeuratObject()

scm_sobj@meta.data <- mrg_tam@meta.data

mrg_tam |> VlnPlot(head(segList$human$human_scSEG, n=3), group.by = 'Platform', pt.size = 0) &
  xlab('Platform')

scm_sobj |> VlnPlot(head(segList$human$human_scSEG, n=3), group.by = 'Platform', pt.size = 0) &
  xlab('Platform')

up_seek <- scm_sobj |>
  FindMarkers(ident.1 = 'TT', ident.2 = 'IT', group.by = 'genotype', logfc.threshold = .1) |>
  as_tibble(rownames = 'gene') |>
  mutate(type = case_when(avg_log2FC > 1 & p_val_adj < .01 ~ 'Up in TT',
                          avg_log2FC < -1 & p_val_adj < .01 ~ 'Down in TT',
                          .default = 'NS'))

up_seek |>
  count(type)

up_seek_name <- up_seek |>
  filter(type != 'NS')
  
up_seek |>
  ggplot(aes(avg_log2FC, -log10(p_val_adj), color = type)) +
  geom_point() +
  geom_vline(xintercept = c(-1,1), linetype = 'dashed') +
  geom_hline(yintercept = 10, linetype = 'dashed') +
  geom_text_repel(data = up_seek_name, aes(label = gene), color = 'black') +
  scale_color_manual(values = c('blue','grey','red'), label = c('Up in IT','NS','Up in TT')) +
  expand_limits(x = c(-3,3)) +
  theme_pubr() +
  labs_pubr() +
  labs(title = 'C1QC+ TAM DEG: TT vs IT')

up_seek <- mrg_tam |>
  FindMarkers(ident.1 = 'TT', ident.2 = 'IT', group.by = 'genotype', logfc.threshold = .1) |>
  as_tibble(rownames = 'gene')

up_seek <- up_seek |>
  mutate(type = case_when(avg_log2FC > .75 & p_val_adj < .01 ~ 'Up in TT',
                          avg_log2FC < -.75 & p_val_adj < .01 ~ 'Down in TT',
                          .default = 'NS'))

up_seek |>
  count(type)

up_seek_name <- up_seek |>
  filter(type != 'NS')

up_seek |>
  ggplot(aes(avg_log2FC, -log10(p_val_adj), color = type)) +
  geom_point() +
  geom_vline(xintercept = c(-.75,.75), linetype = 'dashed') +
  geom_hline(yintercept = 5, linetype = 'dashed') +
  geom_text_repel(data = up_seek_name, aes(label = gene), color = 'black') +
  scale_color_manual(values = c('blue','grey','red'), label = c('Up in IT','NS','Up in TT')) +
  expand_limits(x = c(-3,3)) +
  theme_pubr() +
  labs_pubr() +
  labs(title = 'C1QC+ TAM DEG: TT vs IT')

# process separately -----------
sobj@meta.data %>%
  select(contains(c('Cluster', 'genotype', 'sample'))) ->
  sobj@meta.data

sobj %>%
  FindVariableFeatures() %>%
  ScaleData() %>%
  RunPCA() %>%
  RunUMAP(dims = 1:20) %>%
  FindNeighbors(dims = 1:20) %>%
  FindClusters() ->
  sobj

DimPlot(sobj, group.by = 'seurat_clusters', reduction = 'pca')
DimPlot(sobj, group.by = 'Sub_Cluster')

sg_sobj %>%
  RunUMAP(dims = 1:20) %>%
  FindNeighbors(dims = 1:20) %>%
  FindClusters() ->
  sg_sobj

# identify cell type by public ref ------------
hpca <- read_rds('CRC-I/ref/HumanPrimaryCellAtlas.rds')

monaco <- celldex::MonacoImmuneData()

annotate_cell_type_by_singler <- function(
    srt,
    ref = hpca,
    ref_label = hpca$label.main,
    new_label = 'singler_labels',
    de_method = 'classic'){
  srt %>%
  as.SingleCellExperiment() %>%
  SingleR::SingleR(
    ref = ref,
    labels = ref_label,
    clusters = srt$seurat_clusters,
    de.method = de_method) ->
  singler_sce

  SingleR::plotScoreHeatmap(singler_sce, 
                            show.pruned = TRUE,
                            show_colnames = TRUE)
  
  tibble(seurat_clusters = singler_sce@rownames, singler_labels = singler_sce$pruned.labels) %>%
  rename_with(~new_label, singler_labels) %>%
  right_join(rownames_to_column(srt@meta.data, 'rownames')) %>%
  column_to_rownames('rownames')
}

annotate_cell_type_by_singler(sobj) -> sobj@meta.data

DimPlot(sobj)

annotate_cell_type_by_singler(sg_sobj) -> sg_sobj@meta.data

DimPlot(sg_sobj)
DimPlot(sg_sobj, group.by = 'singler_labels')

# check_duplicates is a latent parameter in RunTSNE, set to FALSE if error is thrown
sg_sobj <- RunTSNE(sg_sobj, check_duplicates = FALSE)

TSNEPlot(sg_sobj, group.by = 'singler_labels')

# remove stromal cells and re-analyze crc0926 ----
sg_sobj <- subset(sg_sobj, singler_labels %in% c(
  'B_cell',
  'DC',
  'Monocyte',
  'NK_cell',
  'T_cells'
))

sg_sobj <- sg_sobj %>%
  FindVariableFeatures() %>%
  ScaleData() %>%
  RunPCA() %>%
  RunUMAP(dims = 1:20) %>%
  FindNeighbors(dims = 1:20) %>%
  FindClusters()

DimPlot(sg_sobj)

sg_sobj@meta.data <- annotate_cell_type_by_singler(
  sg_sobj,
  ref = monaco,
  ref_label = monaco$label.fine,
  new_label = 'monaco_label')

DimPlot(sg_sobj, group.by = 'monaco_label')

sg_sobj@meta.data %>%
  dplyr::count(seurat_clusters, monaco_label)

# prepare for ecotyper -----------
tribble(~seurat_clusters, ~eco_label,
       6 , 'B.cells',
       7 , 'Monocytes.and.Macrophages',
       11 , 'Monocytes.and.Macrophages',
       1 , 'PMNs',
       12 , 'Dendritic.cells',
       8 , 'NK.cells',
       2 , 'PCs',
       9 , 'PCs',
       10 , 'PCs',
       13 , 'PCs',
       14 , 'Dendritic.cells',
       0 , 'CD4.T.cells',
       3 , 'CD8.T.cells',
       4 , 'CD8.T.cells',
       5 , 'CD4.T.cells') -> eco_tbl

sg_sobj@meta.data <- sg_sobj@meta.data %>%
  rownames_to_column('rownames') %>%
  mutate(seurat_clusters = as.numeric(seurat_clusters)) %>%
  left_join(eco_tbl) %>%
  column_to_rownames('rownames')

sg_sobj$genotype <- 'TT'

eco_input <- sg_sobj@meta.data %>%
  rownames_to_column('ID') %>%
  select(c(ID, eco_label, genotype))

colnames(eco_input) <- c('ID', 'CellType', 'Sample')

write_tsv(eco_input, 'CRC-I/results/crc0926_immune_eco_anno.tsv')

sg_sobj@assays$RNA@data %>%
  as.data.frame() %>%
  rownames_to_column('Cell-ID') -> df

write_tsv(df, 'CRC-I/results/crc0926_immune_mat.tsv')

sobj@assays$RNA@data %>%
  as.data.frame() %>%
  rownames_to_column('Cell-ID') %>%
  write_tsv('CRC-I/results/zy2020_tme_mat.tsv')

sobj@meta.data %>%
  rownames_to_column('ID') %>%
  mutate(CellType = case_when(
    str_detect(Sub_Cluster, 'CD4') ~ 'CD4.T.cells',
    str_detect(Sub_Cluster, 'CD8') ~ 'CD8.T.cells',
    str_detect(Sub_Cluster, 'NK') ~ 'NK.cells',
    str_detect(Sub_Cluster, 'Plasma') ~ 'PCs',
    str_detect(Sub_Cluster, 'hB') ~ 'B.cells',
    str_detect(Sub_Cluster, 'DC') ~ 'Dendritic.cells',
    str_detect(Sub_Cluster, 'Mast') ~ 'Mast.cells',
    str_detect(Sub_Cluster, 'hM') ~ 'Monocytes.and.Macrophages',
    TRUE ~ Sub_Cluster
  )) %>%
  select(c(ID, CellType, genotype)) %>%
  rename(Sample = genotype) %>%
  write_tsv('CRC-I/results/zy2020_immune_eco_anno.tsv')
  
# identify cell type by 10x Zhang-Yu 2020 as ref ----------
sobj %>% as.SingleCellExperiment() -> sce10x

new_meta <- annotate_cell_type_by_singler(
  sg_sobj,
  ref = sce10x,
  ref_label = sce10x$Sub_Cluster, 
  new_label = 'Sub_Cluster',
  de_method = 'wilcox')

# save integrated file ---------
write_rds(sg_sobj, 'CRC-I/data/crc0926_singler.rds')
write_rds(sobj, "CRC-I/data/Zhang-Yu-2020/tenx_singler.rds")
write_rds(scM_sobj, 'CRC-I/data/scMerge_CRC221230.rds')